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General relativity predicts the gravitational radiation signatures of mergers of compact binaries, 
such as coalescing binary black hole systems. Derivations of waveform predictions for such systems 
are required for optimal scientific analysis of observational gravitational wave data, and have so far 
been achieved primarily with the aid of the post-Newtonian (PN) approximation. The quality of this 
treatment is unclear, however, for the important late inspiral portion. We derive late-inspiral wave- 
forms via a complementary approach, direct numerical simulation of Einstein’s equations, which has 
recently matured sufficiently for such applications. We compare waveform phasing from simulations 
covering the last ~ 14 cycles of gravitational radiation from an equal-mass binary system of non- 
spinning black holes with the corresponding 3PN and 3.5PN orbital phasing. We find agreement 
consistent with internal error estimates based on either approach at the level of one radian over 
~ 10. cycles. The result suggests that PN waveforms for this system are effective roughly until the 
system reaches its last stable orbit just prior to the final merger. 
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Compact astrophysical binaries spiral together due to 
the emission of gravitational radiation. Calculating the 
dynamics of these systems, and the corresponding gravi- 
tational waveforms, has been a central problem in general 
relativity for several decades. With first-generation inter- 
ferometric detectors such as LIGO, VIRGO, and GEO600 
now operating, and development moving forward on the 
space-based LISA mission, the need for accurate and re- 
liable waveforms for gravitational wave data analysis has 
become urgent. 

Post-Newtonian (PN) methods, based on expansions 
in the parameter e ~ v/c, have been the major analytic 
tool used to calculate the system dynamics and wave- 
forms during the early part of the inspiral, when the 
binary components are relatively widely separated and 
thus have a small orbital frequency [1] . Currently, grav- 
itational wave data analysis for binary inspiral relies on 
waveforms derived from PN methods [2]. The current 
predicted orbital phase is available up to 0(e 7 ), which 
is referred to as 3.5PN order. However the convergence 
properties of the PN sequence are not well understood, 
and it is not yet clear how well PN predictions work late 
in the inspiral when frequencies are high. 

Numerical relativity, in which the full set of Einstein’s 
equations is solved on a computer, is needed to handle 
the final stages of the binary evolution, when the com- 
ponents inspiral rapidly and merge. Recently, there has 
been dramatic progress in the use of numerical relativity 
to simulate the final inspiral and merger of black holes 
[3, 4, 5, 6, 7, 8, 9, 10, 11]. These breakthroughs have al- 
lowed numerical simulations with increasingly wider ini- 
tial separations, producing longer wavetrains. Linking 
such simulations with the PN calculations and comparing 
their waveforms in the late inspiral regime is a pressing 


concern of gravitational wave data analysis. 

We have carried out numerical simulations of a merg- 
ing equal-mass, nonspinning black-hole binary of suffi- 
cient duration to conduct meaningful comparisons with 
the PN approximation. The black holes start on nearly 
circular orbiting trajectories ~ 1200M before merger, 
where M refers to the mass that the system would have 
had when the black holes were still far apart, before ra- 
diative losses were significant. M is related to time by 
M = 5 x 10 -6 s(M/Mq). In this Letter , we quantitatively 
compare crucial phasing information in the our numeri- 
cal simulation waveforms with phasing in PN waveforms, 
finding striking agreement. 

The numerical simulations were performed using the 
moving puncture method [5, 6, 12]. We use fourth-order 
Runge-Kutta time integration, fourth-order accurate fi- 
nite spatial differencing, and second-order-accurate ini- 
tial data. Adaptive mesh refinement is used to resolve 
both the dynamics near the black holes and the propaga- 
tion of the gravitational waves [8]. We performed physi- 
cally equivalent runs at three different maximum resolu- 
tions: low (3M/64), medium (3M/80), and high (M/32). 
We find fourth-order convergence of the Hamiltonian con- 
straint, and better than second-order convergence of the 
momentum constraints during the runs. 

The simulations begin at an angular gravitational wave 
frequency u ~ 0.051M -1 . The frequency then sweeps 
upward through roughly an order of magnitude while 
the black holes undergo ~ 7 orbits, thus producing ~ 14 
gravitational wave cycles before merger. For such long- 
lasting simulations, the primary consideration in provid- 
ing a realistic initial data model is to set up the orbiting 
black holes with minimal eccentricity, as gravitating bi- 
nary systems of comparable-mass objects are expected to 
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FIG. 1: Gravitational strain waveforms from the merger of 
equal-mass Schwarzschild black holes. The solid curve is the 
waveform from the high resolution numerical simulation, and 
the dashed curve is a PN waveform with 3.5PN order phas- 
ing and 2.5PN order amplitude accuracy. Time £ = 0 is the 
moment of peak radiation amplitude in the simulation. 


circularize rapidly through the emission of gravitational 
radiation. We have selected an initial black hole config- 
uration with relatively little eccentricity of less than one 
percent, as measured below. 

Fig. 1 shows a comparison of gravitational wave strain 
generated by our highest-resolution numerical run and 
that predicted by the PN approximation with 3.5PN 
phasing and 2.5PN (beyond leading order) amplitude ac- 
curacy. The waves are based on the dominant l = 2, m — 
2 spin- weighted spherical harmonic of the radiation, and 
represent an observation made on the system’s equato- 
rial plane, where only one polarization component con- 
tributes to the measured strain. The initial phase and 
initial time of the waves have been adjusted so that the 
frequency and phase for each waveform agree early in 
the simulation, at t — — 1000M. We will quantitatively 
study the evident phase agreement below. 

To conduct comparisons with PN calculations, we need 
to extract an instantaneous gauge-invariant polarization 
phase (f> and angular frequency u> from our simulations. 
These are derived from the first time-derivative of the 
gravitational wave strain, which is a robust quantity in 
the numerical data. This frequency corresponds to the 
sweep rate of the polarization angle of the circularly po- 
larized gravitational wave that can be observed on the 
system’s rotation axis. 

We define eccentricity as deviation from an underlying 
smooth, secular trend. We obtain a monotonic “secular” 
frequency-time relation by modeling the frequency u as 
a fourth-order monotonic polynomial o; c (t), plus an ec- 
centric modulation in the waveform angular frequency, 
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FIG. 2: Fit to estimate eccentricity in the high resolution 
run. The dashed (green) curve is the numerically observed 
frequency. The solid (red) curve is a fit to a secular increase 
with time with a sinusoidal perturbation caused by the ec- 
centricity. The dotted (blue) curve is the “de-eccentrized” 
secular trend lj c . 


do>(£), of the form duj{t) = u(t) — c o c (t) — A sin($(£)), 
where $(£) is a quadratic function of time. The fitting 
procedure is illustrated in Fig. 2, which shows the orig- 
inal numerical frequency c j(t), the fit, and uj c (t). Fit- 
ting this data yields A = 8(±1) x 10“ 4 . For Keple- 
rian systems, conserved angular momentum is propor- 
tional to r 2 u) so the (radial) eccentricity corresponds 
to half the fractional amplitude of frequency modula- 
tions: e = A/{2u). In our case the eccentricity starts 
near 0.008, decreasing by a factor of three by the time 
u c M ~ 0.15. We will compare our simulation with 
non-eccentric PN calculations, with the expectation that 
small eccentricities have minimal effect on the important 
underlying secular trend in the rate at which frequency 
sweeps up approaching merger. 

The phasing of the waveform is critical for gravita- 
tional wave observation. For data analysis, the optimal 
methods for both detection and parameter estimation 
rely on matched filtering, which employs a weighted in- 
ner product that can be expressed in Fourier space as 
<h,s>= f(h*{f)s(f) + h(f)s*(f))/S n (f)df, where h is 
the template being used, s is the signal being analyzed, 
and S n is the one-sided power spectral density of the de- 
tector’s noise [13]. A template that maximizes < h,s > 
will provide an optimal filter. While a template with 
time-varying amplitude can emphasize some portions of 
a signal and not others, the crucial factor is the relative 
phasing of the template and signal. The inner product 
will cease to accumulate in sweeping through frequency 
if the template and the signal evolve to be out of phase 
with each other by more than a half-cycle, decreasing the 
effectiveness of the procedure. 

Our key objective is to compare phasing between nu- 
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merical and PN waveforms. To avoid issues with timing 
alignment, we will compare phases as a function of po- 
larization frequency, which corresponds to twice the or- 
bital frequency in the PN case. For circular inspiral this 
frequency should grow monotonically in time, with the 
frequency w c providing a physical reference of the “hard- 
ness” of the tightening binary. 

Circular inspiral phasing information is typically de- 
rived in PN theory by imposing an energy balance re- 
lation to deduce the rate at which cv c evolves from the 
radiation rate at a specified value of uj c [1] . Though not 
strictly derived in the PN context, this physically sensi- 
ble condition currently allows the determination of the 
chirp rate cb c (cu c ) (or something equivalent) up to 3.5PN 
order [14]. From such a relation information about phase 
and time are determined by integrating d(f)/du c = u c /u c 
and dt/du c = 1 /Cj c . The phasing information can be rep- 
resented by any one of several relations between phase, 
frequency and time. Various approaches take the PN- 
expanded representation of one of these relations as the 
PN “result” for waveform phasing [1, 3, 14]. We consider 
both the PN expansion of the chirp rate, numerically in- 
tegrated as needed, as well as the PN expansion of the 
phase. 

For the purpose of comparison with our numerical sim- 
ulations, we invert the monotonic function oj c (t) to obtain 
the phase as a function of frequency: cj)( u c ) = (f>(t(u c )). 
Note that the effect of eccentricity is not removed from 
0, though the “circularized” frequency u c does provide 
the abscissa according to which phases are compared in 
the different treatments. 

Fig. 3 shows the wave phases; here the phases are ad- 
justed by addition of a constant so that each phase van- 
ishes at u c M = 0.054, corresponding to the time 1000M 
before the radiation merger peak in the high resolution 
numerical simulation. The sequence of numerical results 
converges at fourth-order (verified in Fig. 4), allowing us 
to project by Richardson extrapolation to a fifth-order 
accurate result (thick solid line in Fig.3). 

In Fig. 3 we compare the numerical simulation re- 
sults with two versions of PN-phasing, based on either 
the PN expansion of the phase or on the numerically in- 
tegrated PN expansion of the chirp rate. We show each 
of these at 3PN and 3.5PN order. While the 3PN phase 
seems to agree closely with the extrapolated numerical 
phase, the 3.5PN phase moves away, and develops an 
unphysical negative slope after uj c M ~ 0.15. The agree- 
ment of the extrapolated numerical result with the inte- 
grated PN chirp rate improves from 3PN to 3.5PN, with 
the 3.5PN result showing striking agreement up to about 
u c M ~ 0.15. 

We look more quantitatively at the differences among 
the phases in Fig. 4. The solid curve shows the differ- 
ence between our medium and high resolution results, 
while the dotted curve shows the difference between our 
low and medium resolution results, scaled such that for 
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FIG. 3: Gravitational wave phase, in radians, for numeri- 
cal and two versions of PN waveforms. The sequence of nu- 
merical simulation results approaches the solid curve, which 
agrees well with the phase obtained by numerically integrat- 
ing the 3.5PN expansion of the chirp rate u c (wc)- It also 
agrees with the PN expansion of <j){uj c ) at 3PN order though 
the 3.5PN order curve reaches an unphysical maximum value 
near uj c M ~ 0.2. 


fourth-order convergence the curves should superpose. 
This is indeed demonstrated to good approximation. A 
good estimate for the error of the phase in the high resolu- 
tion run is given by its difference from the phase obtained 
by Richardson extrapolation; this comes out to ~ 93% 
of the med-high (solid) curve shown in the figure. Note 
that the cumulative errors in the numerically-generated 
waveforms accrue primarily before the orbital frequency 
t o c M = 0.1. (This makes sense generally since the simu- 
lations spend much more time at lower frequencies.) 

Without monotonic convergence between the 2PN, 
2.5PN and 3PN at the frequencies considered here, it is 
difficult to estimate errors in the PN phase. Nonetheless 
we tentatively take the difference between the integrated 
3PN and 3.5PN chirp rates, shown by the dashed line in 
Fig. 4, as a measure of PN errors. 

The trend in the slope of these error curves indicates 
the rate at which phase error accumulates, as indepen- 
dently estimated within each approach. Fig. 4 sug- 
gests that phasing errors for our high resolution simu- 
lation accumulate more quickly than PN phasing errors 
for u c M < 0.08 (t < — 300M), with the numerical sim- 
ulation phasing being more accurate than PN at higher 
frequencies. In both cases, the phase error accumulates 
to roughly two radians by c o c M ~ 0.2 as the black holes 
begin to plunge together. 

We now address the central objective of this letter, 
a quantitative comparison of numerical and PN phasing 
results. We compare the Richardson-extrapolated phase, 
our best simulation estimate, with the integrated 3.5PN 
chirp rate result. The difference is shown in the solid 
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FIG. 4: Gravitational wave phase error estimates. Differ- 
ences between phasing from the integrated 3.5PN chirp rate 
and Richardson extrapolation from the numerical simulations 
(solid curve) are small, and are consistent with internal error 
estimates for the numerical simulation results and the PN 
sequence. 

curve of Fig. 4. Note that the phase differences be- 
tween the PN and numerical results (solid curve) accu- 
mulate less rapidly than the estimated errors in either 
the high resolution simulation or the 3.5PN results over 
much of the frequency range shown. This indicates that 
the numerical and post-Newtonian predictions appear to 
be converging to a common answer for waveform phase, 
in a frequency regime where the validity of the PN pre- 
dictions could not previously be assessed. The phase dif- 
ference among the best estimates of the numerical and 
PN approaches begins to grow more significantly after 
about l o c M ~ 0.15 (t ~ —70 M) occurring about a cycle 
before the estimated end of the binary’s last stable orbit. 

Through the frequency range 0.054 < u c M < 0.15 
the net phase difference, measured against frequency, 
amounts to less than one radian, a level of error which 
would be tolerable in many gravitational wave data anal- 
ysis applications. For most of the range studied here, the 
frequency-based phase differences between our high res- 
olution simulation and the integrated 3.5PN chirp rate 
result seem larger than the time-based phase differences, 
which are evident in Fig. 1. This is because the sweep 
rate differences have not yet grown to a significant phase 
difference; the difference would become evident if the 
waveform was continued farther. At the higher frequency 
end of the waveform, the time-based differences would 
appear larger, unless the waveform was aligned to agree 
at the end of the waveform. This sensitivity to time- 
alignment makes frequency- based phase comparisons a 
more reliable indicator of phasing differences. 

Our results provide a crucial cross-validation of PN 
waveforms from the late inspiral of binary merger, with 
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results of new long-lasting numerical simulations. These 
simulations have sufficient accuracy to provide a mean- 
ingful comparison with PN waveforms over the last t ~ 
1000M of the coalescence, specifically addressing a bi- 
nary system of equal-mass non-spinning black holes. We 
find phase agreement consistent with internal phase-error 
estimates conducted in each approach, indicating that 
phase accuracies within a few radians are now achievable 
for this part of the coalescence waveform. 

We emphasize, however, that there is still much impor- 
tant work to be done in improving and further assessing 
PN and numerical simulation waveforms. Certainly we 
have only addressed one case in a large parameter space 
of potential binaries, which will inevitably include sys- 
tems, such as rapidly precessing unequal-mass spinning 
binaries, which are harder to treat with present PN and 
numerical techniques. With either approach, even for our 
simple case, a non-negligible amount of phasing error ac- 
cumulates over the range studied, and more will have ac- 
cumulated at lower- frequencies addressable through the 
PN approximation. We expect continuing developments 
in numerical simulations and the pursuit of higher-order 
PN treatments to be crucial for developing a refined un- 
derstanding of coalescence waveforms, which will be cru- 
cial in some data analysis applications for gravitational 
wave observations. 
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